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Abstract 

This is a follow-up of our recently proposed work on pseudopotential calculation (Ref. [21]) 
of atoms and molecules within DFT framework, using cartesian coordinate grid. Detailed results 
are presented to demonstrate the usefulness, applicability of the same for a larger set of species 
(5 atoms; 53 molecules) and exchange-correlation functionals (local, nonlocal). A thorough com- 
parison on total, component, ionization, atomization energies, eigenvalues, potential energy curves 
with available literature data shows excellent agreement. Additionally, HOMO energies for a series 
of molecules show significant improvements by using the Leeuwen-Baerends exchange potential, 
compared to other functionals considered. Comparison with experiments has been made, wherever 
possible. 
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I. INTRODUCTION 



Development of accurate efficient methods for the electronic structure of molecules, solids, 
clusters has been a topic of increasingly intense interest in chemistry, physics, materials sci- 
ence, etc. In the past two decades, density functional theory (DFT) has been applied with 
remarkable success to investigate the structural (such as geometries, vibrational frequen- 
cies, energetics of chemical reactions, etc.) as well as dynamical (such as photoabsorption, 
photoemission, multi-photon ionization, high-order harmonic generation, etc.) character- 
istics of such systems. In essence, the ground-state electronic energy is partitioned as: 
E[p(r)} = T s [p(r)] + U[p(r)] + E xc [p(r)}, where the three terms in the right-hand side denote 
kinetic energy of a system of non-interacting particles of same density p(r), classical electro- 
static energy (including electron-nuclear attraction, electron-electron and nuclear-nuclear 
repulsion), and exchange-correlation (XC) energy respectively. The single-particle wave 

functions are obtained from the Kohn-Sham (KS) equation as: 

1, 



:V 2 + ^ e// ( 



^(r) = e^(r) (1) 
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where v e ff(r) contains the relevant one- and two-body potentials, while the electron density 
is given by the squared sum of the occupied orbitals, p(r) = J2i ./il^M 1 ")! 2 - 

Broadly speaking, two approaches have gainedpopularity for the practical solution of 
above KS equation. In the real-space method [l-5j the discretized KS equation is solved 
iteratively on a mesh using either finite-difference, finite-element or wavelets. The grid-based 
matrix presentation produces structured, highly banded matrices. Moreover, the near locality 
(the potential operator is diagonal in coordinate space whereas the Laplacian operator is 
nearly local) can potentially remove the omnipresent nagging problem arising due to the 
basis-set incompleteness in the contrasting basis-set approaches. Furthermore, they are 
easily amenable to the so-called linear-scaling methods. While such schemes usually require 
a larger number of grids to achieve physically meaningful results, through introduction 
of higher-order and multigrid techniques, the effective grid can be reduced significantly. In 
basis-set approaches, on the other hand, the single-particle wave functions are represented by 
a variety of functions such as Slater-type orbitals, gaussian type functions (GTF), numerical 
functions, plane waves (PW), augmented plane waves, linear muffin-tin orbitals etc. While 
currently, there is a preponderance of the atom-centered localized GTFs in the field of 
molecular quantum chemistry {f| (chiefly due their ability to deal with some of the important 



multi-center integrals analytically), ah initio calculation in the condensed matter regime is 
almost exclusively dominated by the PWs [7| (mainly because the Fast-Fourier transform 
(FFT) technique can be employed to take advantage of the periodicity of such systems). This 
work is concerned with the basis-set approach to tackle the electron-structure calculation of 
atoms and molecules. Recently there has been effort to combine the basis sets; e.g., in jg], a 
gaussian basis set was used to expand the KS molecular orbital (MO), whereas an auxiliary 
augmented PW basis was utilized to express the electronic charge density. 

Excepting a very few attempts such as the grid-free approach {q], that uses a resolution 
of identity to evaluate multi-center integrals over functionals at the expense of an auxil- 
iary basis set, a large majority of basis-set-based molecular DFT calculations employ some 
carefully selected suitable 3D quadrature for a sensible choice of grid points in space. A 
very successful common scheme jlo| involves partitioning the 3D molecular integrand into 
single-center discrete atomic "cells". These can be treated using standard numerical tech- 
niques individually and then summing these up using some appropriate weight functions 
leads to the desired final result. In this so-called atom-centered grid (ACG) method, the 
mono-centric atomic integrals are computed by separating the radial and angular compo- 
nents. The former has typically been performed by introducing several quadratures and 
various mapping schemes. Some of these are Gauss-Chebyshev quadrature of second kind, 
Chebyshev quadratures of first and second kind, Gaussian quadrature, Euler-MacLaurin for- 



mula, numerical quadratures, etc. |llHl5j . Angular integrations are handled usually by the 
Lebedev spherical method 16l-ll8j. In another development, 19] 3D integration of molecular 
integrands were performed numerically based on a division of space and subsequent integra- 
tion over the resulting regions by product Gauss rule. A variational integration mesh 20] , 
which depends on the position of individual atoms has also been reported by breaking the 
space into three different regions, viz., atomic spheres, excluded cubic region and interstitial 
parallelepiped. Other numerical grids [2l] have been proposed as well for gaussian basis set 
calculations. Recently there has been an attempt to connect the cartesian coordinate grid 
(CCG) and ACG by a divided-difference polynomial interpolation which can translate the 
electron density and gradients from the former to the latter {22]. 

In a previous article 23], henceforth referred to as I, efforts were made to employ the 
CCGs in the context of atoms and molecules within the linear combination of GTF ansatz 
of DFT. In this work, respective quantities such as the localized atom-centered basis set, the 
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two-body classical Coulomb repulsion and the non-classical XC potentials were directly set 
up on the 3D CCG. Local-density approximated (LDA) XC functionals of the homogeneous 
electron gas in conjunction with the Hay-Wadt type pseudopotential was used. A multitude 
of quantities like viz., individual energy contributions to total energy, eigenvalues, potential 
energy curves, ionization energies, as well as atomization energies were investigated. The 
relevance and performance of the CCGs were judged on 12 molecules (from small to medium 
size) and 3 atoms. In each of these cases practically identical results as those from the ACG 
and grid-free results were obtained. Now, it is well-known that the local density functionals 
suffer from a number of problems and hence it is essential to use more accurate functionals, 
namely those incorporating the gradient and Laplacian corrections. We have a number 
of objectives in this article. Firstly, we want to assess the validity and efficacy of the 
aforementioned scheme of I for a larger set of species to justify and validate its future usage. 
Secondly its scope of applicability is broadened by incorporating the gradient corrected XC 



functionals (for demonstration, the popular non-local Becke exchange 24j and Lee-Yang- 
Parr (LYP) correlation [25[ is used). First, a detailed test on the convergence of our results 
on 5 atoms and 12 molecules is made as in I for the BLYP XC functional. In the next 
step, additional 41 molecules are considered both at the LDA and BLYP level. Detailed 
comparison with reference theoretical and experimental results are made, wherever possible. 
Besides, we also report the highest-occupied molecular orbital (HOMO) energies obtained 



with the Leeuwen-Baerends (LB) 26|, ]27j] exchange potential for the latter 41 molecules 
(see Section III for motivational details). In all cases, this significantly improves both LDA 
and BLYP HOMO energies. Section II gives a summary of the method of calculation and 
computational aspects. Section III presents a discussion on the results obtained. We end 
with a few concluding remarks about the future prospects in Section IV. 

II. METHODOLOGY AND IMPLEMENTATION 

This section sketches the essential steps involved in the ground-state calculation of a 
many-electron system within the KS DFT framework used in the present work; further 



details can be found in 23 ] 



The KS MOs {V'f ( r )}> cr = a, (3 are linearly expanded in terms of a set of K known basis 
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functions as, 

^) = EW. i = l,2,-~,K (2) 

where the set {Xfj.( r )} denotes the contracted gaussian functions centered on the constituent 
atoms while {Cf^} contains the contraction coefficients for the orbital ipf (r). The individual 
spin-densities are given by, 

p a (r) = e icwr (3) 

where P CT stands for the respective density matrices. 

Substitution of these terms into the energy expression, followed by minimization with 
respect to unknown coefficients C° i7 subject to the orthogonality constraint leads to the 
following matrix KS equation which is akin to the unrestricted Pople-Nesbet equation in 
Hartree-Fock (HF) theory, 

F a C a = SC a e a . (4) 

Here S and F denote the K x K real, square symmetric overlap and total KS matrices 
respectively. C stands for the eigenvector matrix containing the expansion coefficients CZ^ 
while the orbital energies 6{ are embedded in the diagonal matrix e. The KS matrix is 
conveniently written as, 

f; u = H%r + j,* + f* c ° (5) 

where the first two terms in the right-hand side signify the matrices of one-electron bare- 
nucleus Hamiltonian and the classical Coulomb repulsion. The one-electron overlap, kinetic- 
energy and nuclear-attraction integrals involved are the same as encountered _in the GTF- 
based HF schemes and are computed using standard recursion algorithm 
pseudopotential matrix elements in the gaussian basis are taken taken from the GAMESS 
quantum chemistry program output 
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29|. The 



convolution technique 
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301 ]. The Coulomb potential is evaluated using a Fourier 



32j | . which effectively uses a Ewald summation type decomposition 



in terms of short- and long-range contributions, 

1 erf (or) erfc(ar) , „ Wt/ 

- = — i — - H i — - = v l ° ng {r) + vf OIt (r). (6) 

Here erf(x) and erfc(x) identify the error function and complimentary error function respec- 
tively. Note that the latter goes to zero exponentially fast at large r. The parameter a 
controls the range on which <y*; hort is nonzero. After a thorough check on the convergence 



with respect to a, a value within the range of 6-8 seemed satisfactory and we employed 7 for 
all the results reported in this work. This has produced highly satisfyin g re sults for a modest 
number of atomic (3) and molecular (12) systems in our previous work 23j. The short-range 
Fourier integral is calculated analytically, whereas the long-range part is obtained from FFT 
of the real-space values. 

As already mentioned, a major thrust of the current communication is to demonstrate 
the feasibility and viability of our current scheme in the context of so-called "non-local" 
(gradient and Laplacian- dependent) XC functional which would be necessary for future 
chemical applications. For this, we first test this with two of the widely used functionals, 
namely the Becke exchange 24| and LYP correlation 25] (for convenience, an alternative 
equivalent form (33! containing only the first derivative has been mostly used in practice). 
The homogeneous electron-gas correlation of Vosko-Wilk-Nusair (VWN) 34[ is used in all 
the LDA calculations. Following [35(, the gradient- dependent functionals can be treated 
without evaluating the density Hessians by using a finite-orbital basis expansion. To this 
end, XC contributions of the KS matrix is written in the following working form, 



rpXCa 

" [IV 



df_ 



dr 



(7) 



where 7, - 12 - - ^~ - - r "" J 

quantities p a , pp and their gradients. The BLYP functionals are implemented using the 



|Vp a | 2 , 7^/3 = V p a • Vp/?, 7/3/3 = |Vp/3| 2 and / is a function only of the local 



Density Functional Repository program 



36|. The two-electron matrix elements cannot be 



evaluated analytically; here we use direct numerical integrations on the CCG grid. Note 
that some of the existing DFT codes |37J] use an alternative route of fitting these by an 
auxiliary basis set, the so-called discrete variational method 38M40| . The generalized matrix- 
eigenvalue problem is solved using the standard LAPACK routine 4l| accurately efficiently. 
Self-consistent set of MOs, density and energy are obtained in the usual manner subject 
to the convergence of (a) potential (b) total energy and (c) eigenvalues. Tolerance of 10~ 6 
a.u. was used for (b) and (c), while a value of 10~ 5 for (a). Atomic units employed, unless 
otherwise mentioned. 



6 



III. RESULTS AND DISCUSSION 



At first, Table I shows a comparison of our calculated ground-state energy components 
with respect to the number of mesh points N r (r G {x,y,z}) as well as the grid spacing 
h r , for CI2 and HC1 at internuclear distances 4.2 and 2.4 a.u respectively employing the 
BLYP density functional. In this occasion, we adopt a similar presentation strategy as in I. 
A series of calculations were performed in the same fashion to test the convergence of our 
results with respect to the grid parameters. These numerical experiments produced very 
similar conclusions as reached in for LDA XC-case in I. Hence out of eight, we eventually 
report here the results from two such sets only (to avoid too many entries in the table), 
viz., (i) Set A with N r = 64, h r = 0.3, (ii) Set B with N r = 128, h r = 0.2, which suffices 
to illustrate the important points. The reference theoretical results presented in this and 
all other tables throughout the article imply those obtained from the GAMESS quantum 
chemistry program [30]. They use same basis set, effective core potential and employing the 
"grid" option. The corresponding results from the "gridfree" option are quoted in footnotes 
for convenience, in several occasions. Now onwards, we will refer to them as grid and grid-free 
theoretical reference results. The former uses the Euler-McLaurin quadrature for the radial 
integrations and Gauss-Legendre quadrature for the angular integrations. The convergence 
of energies and other quantities with respect to the radial and angular grid was monitored by 
performing two extra set of calculations (i) N r , N e , = 96, 36, 72 (ii) N r , N 9 , = 128, 36, 
72, besides the default grid option (N r , Ng, = 96, 12, 24), where the three integers denote 
the respective number of integration points in r, 6, (ft directions. Note that all the 3 grids 
offered very similar results; for example, out of the 17 atoms and molecules, total energies 
remain same upto 5th decimal place for 8 species for all 3 grids. In the remaining cases, they 
differed slightly among each other as the grid parameters changed; the largest deviation in 
total energy being 0.00064 a.u. for Na2Cl2 and for all others it is well below 0.00007 a.u. 
However, passing from the default grid to (ii) gradually improves N. In this and all other 
tables in the article, we have quoted (ii) results for the reference grid-DFT values. The 
grid-free implementation uses the resolution of identity to simplify the molecular integrals 
enabling their analytical evaluation and obviating the necessity of using grid quadratures. 
Quantities considered are the same as those in I, viz., various energies such as kinetic (T), 
total nucleus-electron attraction (V t ne ), classical Coulomb repulsion (Eh), exchange (E x ), 
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TABLE I: Variation of the energy components and N with respect to the grid parameters for CI2 
and HC1 with reference values. BLYP results in a.u. 







Cl 2 (R = 4.2 a.u.) 






HC1 (R = 2.4 a.u.) 




Set 


A 


B 


Ref. [301 


A 


B 


Ref. [30] 


N r 


64 


128 




64 


128 




h r 


0.3 


0.2 




0.3 


0.2 




(T) 


11.21504 


11.21577 


11.21570 


6.25431 


6.25464 


6.25458 


(vn 


-83.72582 


-83.72695 


-83.72685 


-37.29933 


-37.29987 


-37.29979 


(E h ) 


36.74464 


36.74501 




15.86078 


15.86103 




(E x ) 


-5.29009 


-5.29015 




-3.01023 


-3.01026 




(E c ) 


-0.37884 


-0.37892 




-0.21171 


-0.21174 




(vn 


31.07572 


31.07594 


31.07594 


12.63884 


12.63903 


12.63901 


{Enu) 


11.66667 


11.66667 


11.66667 


2.91667 


2.91667 


2.91667 


(V) 


-40.98344 


-40.98434 


-40.98424 


-21.74382 


-21.74417 


-21.74411 


(Eel) 


-41.43506 


-41.43524 


-41.43522 


-18.40618 


-18.40620 


-18.40620 


(E) 


-29.76840 


-29.76857 


-29.76855 a 


-15.48951 


-15.48953 


-15.48953 6 


N 


14.00006 


14.00000 


13.99998 


8.00002 


8.00000 


8.00000 



°The grid-free DFT value is -29.74755 a.u. 
^The grid-free DFT value is -15.48083 a.u. 



correlation (E c ), total two-electron potential {V t ee ), nuclear repulsion (E nu ), total potential 
(V), electronic (E e i) and total energy (E) respectively, as well as the integrated electron 
density N. Evidently, as in the LDA case, both sets offer excellent agreement in total 
energies and component energies with literature values for both molecules. The individual 
two-body energy terms were not available in the reference output and thus could not be 
directly compared. As expected, for obvious reasons Set B results are closer to reference 
values than Set A, but only marginally. For CI2 this is slightly more pronounced than that 
for HC1. Absolute deviations in total energy for Set B for CI2 and HC1 are only 0.00002 and 
0.00000 a.u. respectively. In both cases, there is slight improvement in N, as we move from 
Set A to B. For all practical purposes, Set A is adequate enough for both of them. Note 
that reference grid-free DFT energies differ substantially from the corresponding grid-DFT 
values. 

To further demonstrate the accuracy and reliability of present results, in Table II, cal- 
culated negative eigenvalues for CI2 and HC1 (using BLYP XC combination) are presented 
at the same geometries of previous table, along with those obtained from reference 30 ]. 
Once again, excellent agreement is observed for both molecules. Sets A and B results match 
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TABLE II: Comparison of the calculated negative eigenvalues of CI2 and HC1 with the reference 
values. BLYP results are given in a.u. 



MO Cl 2 (R = 4.2 a.u.) 

Set A B Ref. [30] 

2cr 9 0.8143 0.8143 0.8143 

2cr u 0.7094 0.7094 0.7094 

3cr 9 0.4170 0.4171 0.4171 

1tt xu 0.3405 0.3405 0.3405 

lTTy U 0.3405 0.3405 0.3405 

lTv xg 0.2778 0.2778 0.2778 

lir vq 0.2778 0.2778 0.2778 



MO HC1 (R = 2.4 a.u.) 



A B Ref. [30] 

2cr 0.7707 0.7707 0.7707 

3cr 0.4168 0.4167 0.4167 

1tt x 0.2786 0.2786 0.2786 

1tt v 0.2786 0.2786 0.2786 



completely with literature values for all the orbital energies except the lone case of 3a g for 
Cl 2 (Set A), where the absolute deviation is only 0.0001 a.u. 

Now Table III tabulates our computed negative total energies of Cl 2 (relative to —29 a.u. 
in columns 2-3, left panel) and HC1 (relative to —15 a.u. in columns 6-7, right panel) for 
Sets A and B using the BLYP XC functionals. A broad range of internuclear distance is 
considered in columns 1 and 5 (3.5-5.0 a.u. for the former and 1.5-3.0 a.u. for the latter) 

n 

and the corresponding grid DFT results obtained from the GAMESS program [30| are listed 
in columns 5 and 8 for comparison. These are depicted vividly for smaller ranges of R in 
Fig. 1. Clearly, for both the molecules, Sets A and B results are practically coincident on 
reference values for the entire range of R. For CI2, maximum absolute deviation is 0.0001 
a.u. with Set B and 0.0002 a.u. (only in 2 instances) with Set A. However, for HC1, the 
two corresponding maximum deviations are 0.0001 a.u., for both sets. This is anticipated 
from the results of Table I, where we noticed results from these two sets confirmed to each 
other more for HC1 than for CI2. This discussion clearly illustrates the faithfulness of current 
calculation. 

For additional test, Table IV reports kinetic, potential and total energies as well as N 
for 15 species (5 atoms and 10 molecules) calculated using the BLYP XC functional. With 
the exception of Mg and S, these are the same species studied in Table V of I, using the 
LDA XC functionals. In this and all other tables henceforth, the experimental geometries 



in the NIST database 



42] are used. The component energies show similar agreements with 



the reference values; hence omitted to avoid crowding. The respective grid-DFT results 
obtained from 30] are presented side by side for comparison. These are ordered in terms of 
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TABLE III: Calculated potential energy curves of CI2 and HC1 for grid Sets A, B, along with 
literature values (grid-DFT). Negative values are given in a.u. 



R (a.u.) CI2 (Total energy relative to —29 a.u.) 



bet 


A 


r> 

L> 


Kci. loUJ 


3.50 


O 

U. 1 Uoz 


U. I Uoz 


U. iOoo 


3.60 


0.7231 


U. (261 


0.7231 


3.70 


O 79QQ 
U. IO0O 


0.7384 


n n*Q O A 
U. / oo4 


3.80 


0.7497 


0.7498 


U. /4y / 


3.90 


0.7579 


0.7580 


0.7579 


a nn 
4.UU 


U. / OOD 


U. /DOO 


U. /DOO 


4.10 


0.7668 


0.7669 


0.7669 


4.20 


0.7684 


0.7686 


0.7685 


4.30 


0.7686 


0.7688 


0.7687 


4.40 


0.7676 


0.7678 


0.7678 


4.50 


0.7658 


0.7659 


0.7659 


4.60 


0.7632 


0.7633 


0.7633 


4.70 


0.7599 


0.7601 


0.7601 


4.80 


0.7565 


0.7566 


0.7566 


4.90 


0.7527 


0.7528 


0.7528 


5.00 


0.7486 


0.7487 


0.7488 



R (a.u.) HC1 (Total energy relative to —15 a.u.) 





A 


B 


Ref. [30] 


1.50 


0.0533 


0.0534 


0.0533 


1.60 


0.1818 


0.1818 


0.1818 


1.70 


0.2767 


0.2767 


0.2767 


1.80 


0.3464 


0.3464 


0.3464 


1.90 


0.3969 


0.3969 


0.3969 


2.00 


0.4328 


0.4328 


0.4329 


2.10 


0.4578 


0.4577 


0.4577 


2.20 


0.4741 


0.4742 


0.4742 


2.30 


0.4842 


0.4842 


0.4842 


2.40 


0.4895 


0.4895 


0.4895 


2.50 


0.4912 


0.4912 


0.4912 


2.60 


0.4903 


0.4903 


0.4903 


2.70 


0.4874 


0.4874 


0.4874 


2.80 


0.4831 


0.4831 


0.4831 


2.90 


0.4778 


0.4778 


0.4778 


3.00 


0.4718 


0.4718 


0.4718 



-0.76 



-0.761 



-0.763 



5? -0.764 



w -0.765 



-0.766 



-0.767 







Set A 


1 






SetB 


1 






Reference 










I 








Jj 

f 






/ 

// 

// 





















4.2 4.3 4.4 4.5 4.6 4.7 
R(a.u.) 



-0.47 



-0.475 



-0.48 



-0.485 



-0.49 




FIG. 1: Potential energy curves for CI2 (left panel) and HC1 (right panel) for grid Sets A, B. 
Reference grid-DFT results are also given for comparison. 



increasing N as descending the table. First 10 of these use the same grid parameters as in 
Table V of I, i.e., A^ r = 64, h r = 0.4, whereas for last 5 we use iV r = 128, h r = 0.3. Overall, 
the agreement of our results with the reference is excellent. In several occasions (such as 
Na 2 , P, As, Na 2 Cl 2 ), the total energies completely match with them. The largest absolute 
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TABLE IV: Comparison of the kinetic ((T)), potential ((V")), total {E) energies and N for several 
atoms and molecules with the reference grid-DFT results [30j]. BLYP results in a.u. PW=Present 
Work. See text for details. 



System (T) -(V) -(E) N 





PW 


Ref. [301 


PW 


Ref. [301 


PW 


Ref. [30] 


PW 


Ref. [30] 


Na 2 


0.14723 


0.14723 


0.52871 


0.52871 


0.38148 


0.38148 


1.99999 


2.00000 


Mg 


0.24935 


0.24935 


1.06017 


1.06017 


0.81082 


0.81083 


1.99999 


1.99999 


NaH 


0.60093 


0.60088 


1.34698 


1.34693 


0.74606 


0.74605 


1.99997 


1.99999 


P 


2.38891 


2.38890 


8.78249 


8.78248 


6.39358 


6.39358 


4.99999 


4.99999 


As 


2.10490 


2.10494 


8.14207 


8.14211 


6.03717 


6.03717 


4.99999 


4.99999 


S 


3.73143 


3.73145 


13.73598 


13.73599 


10.00455 


10.00454 


6.00000 


5.99999 


Br 


4.27022 


4.27043 


17.36122 


17.36148 


13.09100 


13.09105 


7.00000 


6.99999 


NaCl 


5.83959 


5.83957 


21.01698 


21.01694 


15.17739 


15.17737 


8.00003 


7.99999 


H 2 S 


4.98071 


4.98066 


16.21919 


16.21913 


11.23848 


11.23846 


8.00000 


7.99999 


PH 3 


4.18229 


4.18224 


12.40101 


12.40096 


8.21871 


8.21872 


7.99999 


7.99999 


Br 2 


8.66533 


8.66527 


34.88603 


34.88596 


26.22070 


26.22069 


13.99999 


14.00000 


H2S2 


8.88238 


8.88240 


30.16535 


30.16538 


21.28297 


21.28298 


13.99999 


13.99999 


MgCl 2 


11.75947 


11.75999 


42.54049 


42.54103 


30.78102 


30.78104 


16.00004 


15.99999 


Na 2 Cl 2 


11.68815 


11.68842 


42.12410 


42.12438 


30.43595 


30.43595 


16.00002 


15.99999 


SiH 2 Cl 2 


14.14948 


14.14945 


49.04463 


49.04461 


34.89515 


34.89516 


19.99999 


20.00000 



deviation in total energy is only 0.0013% (for NaH). 

Table V displays the calculated negative HOMO energies, — enoMO ( m a - u -) an d estimated 
atomization energies (in kcals/mole) for the same 12 molecules in Table V of I, using the 
BLYP XC density functional. These are compared with the theoretical [30j as well the 
experimental results from the NIST database 43J. Same grid parameters as those in the 
previous table are employed for these. The — enoMO values completely agree with those 
from the theoretical literature results |30|], except in two occasions (MgCl2, N^C^), where 
the absolute discrepancy remains only 0.0001 a.u. The calculated atomization energies also 



show very good agreement with the reference theoretical values [30j . For Na2 HC1, NaCl and 
H 2 S, our values are identical as those from reference. The largest discrepancy (0.248%) is 
observed for Br2. However, both quantities differ significantly from the experimental values; 
but that is a separate issue (not directly related to the main theme of this work). These 
could possibly be improved further by employing more appropriate basis functions and/or 
better XC functionals, and may be considered in future works. 

The previous work of 23j as well as the ongoing discussion amply demonstrates that the 
approach can produce good-quality results for both LDA and gradient-corrected non-local 
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TABLE V: Comparison of — eHOMo( a - u -) an d atomization energies (kcals/mole) for molecules with 
the literature data using BLYP XC functional. PW=Present Work. 1 a.u.= 627.509471 kcals/mole. 
See text for details. 



System 




— e H OMO (a.U.J 






Atomization energy (kcals/ 


mol) 


PW 


theory |30| 


Expt. [43] 


PW 


theory |30) 


Expt. [13] 


Na 2 


0.1002 


0.1002 


0.1798 


11.51 


11.51 


17.0 


NaH 


0.1421 


0.1421 




45.37 


45.36 


47.2 


HC1 


0.2785 


0.2785 


0.4683 


89.88 


89.88 


102.2 


NaCl 


0.1733 


0.1733 


0.3381 


88.75 


88.75 


97.4 


H 2 S 


0.2190 


0.2190 


0.3843 


156.59 


156.59 


173.2 


PH 3 


0.2287 


0.2287 


0.3627 


218.72 


218.73 


227.1 


Cl 2 


0.2652 


0.2652 


0.4219 


22.89 


22.90 


57.2 


Br 2 


0.2451 


0.2451 


0.3865 


24.28 


24.22 


45.4 


H 2 S 2 


0.2288 


0.2288 


0.3418 


181.66 


181.68 


229.6 


MgCl 2 


0.2697 


0.2698 




164.04 


164.07 


187.4 


Na 2 Cl 2 


0.2020 


0.2021 




228.43 


228.46 


243.1 


SiH 2 Cl 2 


0.2836 


0.2836 


0.4300 


287.92 


287.95 


341.8 



XC functionals for atoms and small-to-medium molecules, within the specifics of basis set. 
This is further validated in Table VI where both LDA and BLYP results are reported for 41 
extra molecules to illustrate the broad sco pe its applicability. The experimental geometries 
are again taken from the NIST database |42|. The kinetic, potential, total energy and N 
are presented. Since it has been clearly established that our obtained results are sufficiently 
accurate, here we omit the reference theoretical values. Only some random checks were 
made to ensure that this is indeed the case. First 29 molecules from top were treated using 
N r = 64, h r = 0.4, while for the remaining 12, we used N r = 128, h r = 0.3 grid. On the light 
of all the results and discussion so far, we believe that these represent the correct values. 

Finally Table VII compares the calculated — €homo (in a.u.) and atomization energies 
(in kcals/mole) with the literature data, for the same molecules at same geometries as in 



previous table. Experimental results, wherever available, are quoted from [43j. In this case 
also, reference results are omitted to avoid crowding. But as expected, they show very little 
deviation from ours. An asterisk in the experimental atomization energies denote the values 
at 298°K; otherwise they imply 0°K values. Both LDA and BLYP results are given for 
these quantities side by side. However for the former we also include results obtained with 
another exchange functional for the following reason. It is well-known that the XC poten- 
tials derived from the simplest LDA or generalized gradient approximation (GGA) suffer 
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TABLE VI: Kinetic ((T)), potential ((V)), total (E) energies and N for several molecules. LDA 
and BLYP results are given in a.u. PW=Present Work. See text for details. 



System (T) -(V) -(E) N 
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SiCl 2 


12.71813 


12.87162 


46.37456 


46.57075 


33.65644 


33.69913 


17.99999 


17.99999 


S3 


11.45897 


11.59948 


41.56900 


41.72809 


30.11003 


30.12861 


17.99999 


17.99999 


CIS2 


13.29292 


13.45889 


48.23575 


48.43954 


34.94284 


34.98065 


19.00000 


19.00000 


P 4 


10.04128 


10.15505 


35.88401 


35.95566 


25.84274 


25.80061 


19.99999 


19.99999 


AICI3 


17.72370 


17.94164 


64.51025 


64.79488 


46.78655 


46.85324 


23.99999 


23.99999 


S2CI2 


19.02565 


19.26416 


68.83640 


69.12631 


49.81075 


49.86216 


25.99999 


25.99999 


PCI3 


19.50036 


19.74348 


70.61349 


70.91249 


51.11314 


51.16902 


26.00000 


26.00000 


S1HCI3 


19.09265 


19.34355 


68.26422 


68.57817 


49.17157 


49.23462 


25.99999 


25.99999 


SiCU 


24.19691 


24.50296 


87.70251 


88.07579 


63.50560 


63.57282 


32.00000 


32.00000 


PCI5 


30.92116 


31.33096 


111.74539 


112.22809 


80.82423 


80.89713 


40.00000 


40.00000 
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from improper asymptotic long-range behavior. Consequently, whereas the ground-state 
total energies of atoms, molecules, solids can be predicted quite satisfactorily by using these 
functionals, the ionization energies obtained via the HOMO energies (usually off by 30-50% 
of the experimental values) as well as the excited states are described rather poorly. As 



mentioned in 



23], our primary objective is to extend this scheme towards the dynamical 



studies of atoms and molecules under the influence of strong field such as a laser (through 
multi-photon ionization, high-order harmonic generation and other related phenomena) via 
TDDFT, that can potentially exploit the remarkable developments made in basis-set DFT 
through many pioneering works over the years. It is a necessary prerequisite that both the 
ionization potential and higher levels be described more accurately. Recently, the modified 
Leeuwen-Baerends (LB) potential 2t]|' v^ a (a,/3 : r), containing two empirical parame- 
ters have been shown to be quite successful in dealing with the above dynamical situations 
of atoms, molecules (see, for example, [4J], and the references therein) as well as for the 
static property calculations including TDDFT-based excited states of molecules. This is 
conveniently written as, 

v LBa (a B-r)= av LDA (r) + v LDA (r) + /^( r )Pa /3 ( r ) , g) 

where a signifies up, down spins and the last term containing gradient correction is 
reminiscent of the popular Becke exchange energy density functional 24J, x a {r) = 
|Vp -(r)|[p .(r)]~ 4//3 is a dimensionless quantity, a = 1.19, ft = 0.01. This ensures the de- 
sired long-range property, i.e., v^ a (r) — > — 1/r, r — > 00. The HOMO energies, obtained 
from LBVWN (LB exchange and VWN correlation) combination are presented in column 
4. It is abundantly clear that LBVWN results are significantly improved over either the 
LDA or BLYP cases. The LDA ionization energies are lower than BLYP values for all 
the molecules considered and LBVWN values are substantially lower than both these two. 
Evidently in our future work on TDDFT as mentioned above, this feature of LB potential 
will be highly exploited. Now columns 6, 7, 8 give the computed LDA, BLYP atomiza- 
tion energies and their experimental analogs. Here also both LDA and BLYP results show 
considerable deviation from the experimental values, which include zero-point vibrational 
corrections as well as relativistic effects. In many cases, LDA results are apparently better 
than their BLYP counterparts. However this observations should not be misconstrued to 
lead to the conclusion that former is a better candidate than the latter. We note that the 
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current work employs the Hay-Wadt-type ab-initio effective core potentials which are more 
suitable for the HF-type approaches. Among other factors, use of pseudopotentials and 
basis sets which are more appropriate for the DFT approaches, may alleviate some of the 
discrepancies encountered here and an undertaking along this direction may be initiated in 
future. Furthermore, large deviations in atomization energies have also been found in other 
recent DFT works involving all-electron calculations and more extended basis sets as well 



(see for example 



451]) . In any case, this is an evolving process and does not interfere with 



the main objective of the present work. Also note that there may be some cancellation of 
errors in the LDA case. Finally note that the extension of this approach to the "all-electron" 
atomic, molecular calculations as well as for very large systems would be relatively difficult 
compared to the present pseudopotential case in terms of the grid requirement, because of 
the presence of extra core electrons. Nevertheless reasonable full calculation of small to 
medium molecules are possible. This is suggested from some of our preliminary studies in 
this direction which I am currently engaged into and may be considered in some future 
communication. 



IV. CONCLUDING REMARKS 



Pseudopotential density functional calculations were performed for atoms and molecules 
within the LCGTF framework using CCG in conjunction with an accurate, efficient Fourier 
convolution technique to represent the classical Hartree potential in real grid. In essence, 
our previous work (I) on the LDA XC functionals has been extended to test its performance 
and validity in the case of gradient-corrected XC functionals which would be necessary for its 
further applications to realistic physical situations. For this purpose, the widely used BLYP 
XC functional was chosen. The calculated results of a variety of quantities such as energy 
components, eigenvalues, potential energy curves, ionization energies, atomization energies 
clearly reveal that, they are practically of the same quality as obtained from the available 
theoretical methods. Furthermore, companion LDA and BLYP calculations were performed 
for a large number of molecules (41) to illustrate its scope of applicability for a broad range 
of systems. Comparison with experiments has been made wherever possible. Additionally, 
for all the molecules studied, the LBVWN results show significant improvements in the 
HOMO energies. This has important relevance to our prospective future works on studying 
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TABLE VII: Negative HOMO energies, — enoMO (i n a.u.) and atomization energies (kcals/mol) for 



more molecules. LDA, LBVWN (LB+VWN), BLYP results are compared with experiment 



43J. 



An asterisk indicates 


298°K values. 


Otherwise 


0°K results 


are given. 


See text for details. 


Molecule 




-«HOMO (a.u.) 




Atomization energy (kcals/mol) 


LDA 


BLYP 


LBVWN 


Expt. [43] 


LDA 


BLYP 


Expt. [43] 


Mg 2 


0.1563 


0.1530 


0.2316 




4.04 


0.22 


0.9 


MgH 2 


0.2248 


0.2221 


0.3471 




113.75 


109.09 




A1H 


0.1741 


0.1715 


0.2836 




72.65 


68.31 


70.3 


SiH 


0.1647 


0.1597 


0.2849 


0.2900 


74.47 


69.59 


70.4 


A1H 2 


0.1655 


0.1631 


0.2755 




127.69 


118.90 


114.7 


Al 2 


0.1407 


0.1400 


0.2371 


0.1984 


22.92 


21.42 


37.0 


PH 


0.2214 


0.2133 


0.3519 


0.3730 


72.05 


67.14 


66.3 


SiH 2 


0.2056 


0.2027 


0.3340 


0.3278 


155.85 


143.93 


144.1 


SH 


0.2229 


0.2174 


0.3736 


0.3830 


84.83 


74.85 


83.9 


HSe 


0.2117 


0.2057 


0.3534 


0.3618 


79.44 


70.16 




PH 2 


0.2170 


0.2111 


0.3504 


0.3610 


152.77 


139.92 


149.2 


SiH 3 


0.2017 


0.1969 


0.3278 


0.2990 


193.37 


171.26 


212.2 


HBr 


0.2688 


0.2603 


0.4147 


0.4292 


91.60 


79.11 


87.5* 


MgS 


0.1850 


0.1766 


0.2882 




55.90 


42.15 


71.7 


NaBr 


0.1818 


0.1729 


0.3057 


0.3050 


87.47 


78.94 


86.8* 


KC1 


0.1481 


0.1419 


0.2805 


0.3859 


96.61 


88.02 


80.3* 


KBr 


0.1528 


0.1449 


0.2760 


0.2903 


87.39 


79.35 


69.8* 


H 2 Se 


0.2152 


0.2075 


0.3524 


0.3635 


167.04 


146.81 




HI 


0.2518 


0.2432 


0.3824 


0.3817 


82.82 


72.07 


45.8* 


SiH 4 


0.3188 


0.3156 


0.4624 


0.4042 


339.43 


312.02 


302.6 


A1S 


0.2342 


0.2240 


0.3497 


0.3491 


93.74 


77.07 


88.4 


MgCl 


0.1848 


0.1804 


0.2750 


0.2753 


78.76 


68.47 


76.5 


P 2 


0.2590 


0.2526 


0.3877 


0.3870 


96.18 


81.73 


116.0 


SiS 


0.2565 


0.2507 


0.3828 


0.3870 


134.04 


114.54 


147.2 


A1C1 


0.2255 


0.2204 


0.3383 


0.3454 


115.69 


100.65 


119.2 


PS 


0.1695 


0.1645 


0.3054 


0.3307 


81.23 


63.43 


105.2 


S 2 


0.2007 


0.2023 


0.3443 


0.3438 


56.75 


52.47 


100.8 


Se 2 


0.1952 


0.1951 


0.3283 


0.3160 


58.40 


54.78 




PCI 


0.2093 


0.2023 


0.3490 




65.12 


49.37 


71.7 


BrCl 


0.2623 


0.2537 


0.4133 


0.4079 


44.95 


25.41 


51.5 


SiH 3 Cl 


0.2780 


0.2704 


0.4317 


0.4267 


337.97 


300.07 


321.7 




0.2514 


0.2448 


0.3909 


0.3804 


190.40 


155.11 


202.7 


s 3 


0.2392 


0.2294 


0.3805 




116.84 


72.14 




C1S 2 


0.2225 


0.2161 


0.3712 




115.08 


73.52 


141.0 


p 4 


0.2712 


0.2575 


0.3964 


0.3432 


200.77 


142.99 


285.9 


AIC13 


0.3081 


0.2976 


0.4603 


0.4414 


278.02 


232.88 


303.4 


S 2 C1 2 


0.2603 


0.2499 


0.4107 


0.3550 


151.27 


90.54 


192.2 


PC13 


0.2747 


0.2660 


0.4266 


0.3638 


189.35 


133.20 


229.5 


SiHCl 3 


0.3076 


0.2971 


0.4632 




335.99 


273.66 


361.3 


S1CI4 


0.3194 


0.3085 


0.4758 


0.4333 


333.91 


258.60 


378.6 


PC1 5 


0.2825 


0.2722 


0.4397 


0.3748 


246.22 


145.33 


303.2 



T6- 



real-time dynamics of atoms and/or molecules in a strong laser field. Incorporation of other 
pseudopotentials more suited to DFT as well as more extended and elaborate basis sets 
would be among some of the important issues which may be considered in recent future. 
More accurate XC functionals could also be employed depending upon the physical system 
concerned and the nature of the problem dealt with. Applications to weakly bonded systems, 
clusters and of course, to larger systems would further consolidate its success. Finally 
although one could think of some inherent errors associated with the incompleteness of the 
grid, this study confirms that with a judicious choice of the grid coupled with a correct 
treatment of the Coulomb potential, these can be reduced to tolerable minima. Thus very 
satisfactory results could be obtained. 

Acknowledgments 

The project was initiated at Prof. D. Neuhauser's laboratory at the Univ. of Califor- 
nia, Los Angeles, where the core of the program was written. It was further extended at 
Prof. S. I. Chu's laboratory at the University of Kansas. I thank them for stimulating 
discussions and providing the computational facilities. Numerous useful discussions with 
Dr. E.I. Proynov is gratefully acknowledged. Warm hospitality offered by UCLA and the 
Univ. of Kansas is greatly appreciated. An anonymous referee is thanked for valuable 
constructive comments. 



[1] S. R. White, J. W. Wilkins and M. P. Teter, Phys. Rev. B 39, 5819 (1989). 

[2] J. R. Chelikowsky, N. Troullier, K. Wu and Y. Saad, Phys. Rev. B 50, 11355 (1994). 

[3] E. L. Briggs, D. J. Sullivan and J. Bernholc, Phys. Rev. B 52, R5471 (1995). 

[4] F. Gygi and G. Galli, Phys. Rev B 52, R2229 (1995). 

[5] T. L. Beck, Rev. Mod. Phys. 72, 1041 (2000). 

[6] T. Helgaker, P. Jorgensen and J. Olsen, Molecular- Electronic Structure Theory (John- Wiley 
& Sons Ltd., 2000). 

[7] M. C. Payne, M. P. Teter, D. C. Allan, T. A. Arias and J. D. Joannopoulos, 
Rev. Mod. Phys. 64, 1045 (1992). 



17 



[8] M. Krack and M. Parinello, Phys. Chem. Chem. Phys. 2, 2105 (2000). 

[9] Y. C. Zheng and J. Almlof, Chem. Phys. Lett. 214, 397 (1993). 

[10] A. D. Becke, J. Chem. Phys. 88, 2547 (1988). 

[11] O. Treutler and R. Ahlrichs, J. Chem. Phys. 102, 346 (1995). 

[12] M. M. Mura and P. J. Knowles, J. Chem. Phys. 104, 9848 (1996). 

[13] C. W. Murray, N. C. Handy and G. J. Laming, Mol. Phys. 78, 997 (1993). 

[14] P. M. W. Gill, B. G. Johnson and J. A. Pople, Chem. Phys. Lett. 209, 506 (1993). 

[15] R. Lindh, P. -A. Malmqvist and L. Gagliardi, Theor. Chem. Ace. 106, 178 (2001). 

[16] V. I. Lebedev and Zh. Vychisl, Mat. Mat. Fiz. 15, 48 (1975). 

[17] V. I. Lebedev and Zh. Vychisl, Mat. Mat. Fiz. 16, 293 (1976). 

[18] V. I. Lebedev and A. L. Skorokhodov, Russ. Acad. Sci. Dokl. Math. 45, 587 (1992). 

[19] P. M. Boerrigter, G. Te. Velde and E. J. Baerends, Int. J. Quant. Chem. 33, 87 (1988). 

[20] M. R. Pederson and K. A. Jackson, Phys. Rev. B 41, 7453 (1990). 

[21] X. Chen, J.-M. Langlois and W. A. Goddard III, Phys. Rev. B 52, 2348 (1995). 

[22] J. Kong, S. T. Brown and L. Fiisti-Molnar, J. Chem. Phys. 124, 094109 (2006). 

[23] A. K. Roy, Int. J. Quant. Chem. 108, 837 (2008). 

[24] A. D. Becke, Phys. Rev. A 38, 3098 (1988). 

[25] C. Lee, W. Yang and R. G. Parr, Phys. Rev. B 37, 785 (1988). 

[26] R. van Leeuwen and E. J. Baerends, Phys. Rev. A 49, 2421 (1994). 

[27] P. R. T. Schipper, O. V. Gritsenko, S. J. A. van Gisbergen and E. J. Baerends, 
J. Chem. Phys. 112, 1344 (20001). 

[28] S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986). 

[29] L. E. McMurchie and E. R. Davidson, J. Comput. Phys. 26, 218 (1978). 

[30] M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H. Hensen, 
S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis and J. A. Mont- 
gomery, J. Comput. Chem. 14, 1347 (1993). 

[31] G. C. Martyna and M. E. Tuckerman, J. Chem. Phys. 110, 2810 (1999). 

[32] P. Minary, M. E. Tuckerman, K. A. Pihakari and G. J. Martyna, J. Chem. Phys. 116, 5351 
(2002). 

[33] B. Miehlich, A. Savin, H. Stoll and H. Preuss, Chem. Phys. Lett. 157, 200 (1989). 

[34] S. H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58, 1200 (1980). 

18 



[35] J. A. Pople, P. M. W. Gill and B. G. Johnson, Chem. Phys. Lett. 199, 557 (1992). 

[36] Density Functional Repository, Quantum Chemistry Group, CCLRC Daresbury Laboratory, 

Daresbury, Cheshire, WA4 4AD, UK. 
[37] J. Andzelm and E. J. Wimmer, J. Chem. Phys. 96, 1280 (1992). 
[38] E. J. Baerends, D. E. Ellis and P. Ros, Chem. Phys. 2, 41 (1973). 
[39] H. Sambe and R. H. Felton, J. Chem. Phys. 62, 1122 (1975). 

[40] B. I. Dunlap, J. W. D. Connolly and J. R. Savin, J. Chem. Phys. 71, 4993 (1979). 

[41] E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Green- 

baum, S. Hammarling, A. McKenney, D. Sorensen, LAPACK Users' Guide. 3rd Ed. (SIAM, 

Philadelphia, 1999). 

[42] R. D. Johnson III (Ed.), NIST Computational Chemistry Comparison and Benchmark 
Database, NIST standard Reference Database. Number 101. Release 14, Sept. 2006. 
(http://srdata.nist.gov/cccbdb ). 

[43] H. Y. Afeefy, J. E. Liebman and S. E. Stein, "Neutral Thermochemical Data" in NIST Chem- 
istry WebBook, NIST Standard Reference Database Number 69, Eds. P. J. Linstrom and 
W. G. Mallard, June 2005, National Institute of Standards and Technology, Gaithersburg 
MD, 20899. flhttp://webbook.nist.gov] ). 

[44] S. I. Chu, J. Chem. Phys. 123, 062207 (2005). 

[45] M. Cafiero, Chem. Phys. Lett. 418, 126 (2006). 



19 



